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We calculate the lattice dispersion relation for three dimensional simulations of scalar fields. We 
argue that the mode frequency of scalar fields on the lattice should not be treated as a function of 
the magnitude of its wavevector but rather of its wavevector decomposition in Fourier space. Fur- 
thermore, we calculate how the lattice dispersion relation differs depending on the way that spatial 
derivatives are discretized when using finite difference methods in configuration space. For applica- 
tions that require the mode frequency as an average function of the magnitude of the wavevector, 
we show how to calculate the radially averaged lattice dispersion relation. Finally, we use the pub- 
licly available framework LATTICEEASY to show that wrong treatment of dispersion relations in 
simulations of preheating leads to an inaccurate description of parametric resonance, which results 
in incorrect calculations of particle number densities during thermalization after inflation. 

o 

I. INTRODUCTION 

^sO The application of numerical simulations in theoretical physics has seen a tremendous surge in the last decade. 

Problems that appeared intractable a few years ago have now been analyzed in depth through advances in hardware 
design as well as extensive research in numerical methods. Computer programs have been employed to solve problems 
ranging from amplitude calculations in quantum field theory [TJ [5] to black hole simulations in numerical general 
relativity [3]. Particularly when there are nonlinear processes involved, analytical approximations often fail to capture 
O | all the feautures of the underlying theories and the use of numerical methods becomes imperative. 

The theory of reheating after inflation [4] has been comprehensively studied using numerical tools which underlined 
the significance of non-perturbative effects towards thermalization [SJ [6| . More recently, considerable interest has been 
directed towards the study of long-lived, coherent objects called oscillons [7HPJ and their relevance in cosmological 

t-H settings [TDTIT2] with recent focus in the period of reheating [T5Ul5| . In this case, analytical techniques prove to be 
unsuccessful to account for the spontaneous emergence of oscillons and numerical simulations provide us with the 
most effective tool to investigate their dynamics. These studies are generally centered around the evolution of scalar 
fields, but models including vector fields have been employed to examine similar processes in Abelian and non-Abelian 
Higgs models (TBI H7j . The interest in such simulations of scalar fields has resulted in the release of publicly available 
numerical frameworks, most notably LATTICEEASY [18] and DEFROST [19] . as well as the more recent PSpectRe 



In this study, we examine how three dimensional dispersion relations are treated in computer programs that use 
finite difference methods, including LATTICEEASY and DEFROST. Even though the aforementioned programs deal 
• • with discretized space, the dispersion relations employed are not adjusted for finite grid size effects. Meanwhile, they 
are important in the setup and evolution of the scalar fields, as they control the initial amplitude of fluctuations in the 
beginning of the simulation, as well as in the definition of the particle number density. As such, we investigate how 
results of preheating simulations are altered using the correct lattice dispersion relations, focusing on previous work 
that is based on LATTICEEASY [5]. Because LATTICEEASY and DEFROST utilize different methods to discretize 
the equations of motion, the lattice dispersion relation in each case differs. We show how to correctly calculate them 
in both approaches. We finally note that even though we focus on preheating simulations, dispersion relations have 
to be adapted to account for lattice effects in any scalar field simulation using finite different methods [2lTl23j . 

This paper is organized as follows: In Section [TT] we derive the lattice dispersion relation for a free massive scalar 
field using the discretization scheme used in LATTICEEASY. In Section |III| we repeat our calculation for a more 



general, isotropic scheme which is used in DEFROST and compare it to LATTICEEASY. In Section IV we show how 
to calculate a radially averaged lattice dispersion relation which is useful whenever mode frequencies are required as 
a function of the wavevector magnitude. Section [V] contains the results that we get from simulations of preheating 
in a X(j) 4 chaotic model of inflation using LATTICEEASY with both the continuous and lattice adjusted dispersion 



relation. We conclude with a summary of our work in Section VI 
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II. LATTICE DISPERSION RELATION USING AN ANISOTROPIC DISCRETIZATION STENCIL 



We consider a free massive scalar field 0(x, t) of mass m in (3+f )-dimcnsional flat Minkowski spacetime with 
Lagrangian density 



2 v /i-r, 2 

Using h — c—1, the equation of motion satisfied by <j) is 



C = hd^f-\m 2 <j> 2 . (1) 



4> - V 2 + m 2 4> = 0, (2) 

where an overdot denotes derivative with respect to time. We uniformly discretize space in a lattice with N 3 points, 
lattice spacing Ax = Ay — Az = a and employ a second-order accurate expression for the Laplacian to get the 
discretized equation of motion for the field at site (i, j, k) 



7 <j>(i+l)jk + <f>(i-l)jk — tyijk 4>i(J+l)k + <t>iU-l)k ~ tyijk 4>ij{k+l) + <t>ij(h-\) — 2<t>ijk 2 , n /o^ 

<Pijk 9 9 9 r m <pijk = U. (6) 

a z a 1 a A 

We note that this discretization of the Laplacian uses the six nearest neighbors of a point which can lead to anisotropic 
propagation of errors in the field evolution. This is the discretization scheme used in LATTICEEASY. We seek a 
solution to Eq. [3] of the form 



<fi(xi jk ,t) = Ae^*^-^, (4) 

where x^ = (ia,ja,ka) and k = (k x ,k v ,k z ). We impose periodic boundary conditions which restricts the allowed 
values of k to k = (k x , k y , k z ) = ^(n x ,n y , n z ) where rii are integers n,; = —N/2 + 1 . . . N/2 and L — Na. Plugging 
Eq. [3] into Eq. [3] and cancelling common factors we get, after some algebra, 



J - A ( sin 2 + Sin 2 hi + sin 2 hi \ _ m 2 = (5) 

a 2 V 2 2 2 ' 

which gives the lattice dispersion relation for this discretization scheme 



to 2 (k x , k y , kg) = feg ff + m 2 



(6) 



with 



e - h 

a z 



2 ® 



(7) 



Unlike the continuous dispersion relation a; 2 = k 2 + m 2 , the lattice dispersion relation formulates the mode frequency 
w as a function of the vector decomposition of k and not the magnitude k = |k|. Two modes ki = (k x ,k y ,k z ) and 
k2 = {k*, ky,k* z ) can have |ki| = |k 2 | but w(ki) ^ w(k 2 ). For small values of ki, Eq. [olreduces to the continuous limit 
which is expected since large wavelength modes are not significantly affected by finite grid size effects. Large values 
of ki, however, give mode frequencies that are much different to what one would get using uj 2 = k 2 + m 2 . 

The validity of Eq. [6] is demonstrated by solving Eq. [3] on a 3d lattice. We use N = 128, a — 0.5, m = 1 and 
propagate </>(x, t) in time using a symplectic second-order Velocity- Verlet algorithm. We initialize the field in Fourier 
space with random amplitudes for each mode and then track the evolution of two modes with the same wavevector 
magnitude. We have chosen the modes with ki = ^(40, 28, 25) and k 2 = ^(53,14,2) for illustration. They both 
have | kx | = |k 2 | = 2 ^v / 3009. The results are summarized in Table I Even though both wavevectors have the same 
magnitude, their numerically calculated frequencies are significantly different. The continuous dispersion relation 
fails in the numerical analysis, assigning to both modes the same incorrect frequency. The data also shows excellent 
agreement between Eq. [6] and the numerically calculated frequencies. 
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Wavevector 




u) — \/k 2 + m 2 


Numerical 


£(40,28,25) 


4.8791m 


5.4774m 


4.8707m 


£(53,14,2) 


4.2091m 


5.4774m 


4.2169m 



TABLE I: Mode frequencies for two modes computed numerically, their values in the continuous case cj 2 = k 2 + m 2 and when 
calculated using Eq. [5] 

III. LATTICE DISPERSION RELATION USING ISOTROPIC DISCRETIZATION 

The Laplacian term in Eq. [2] can be discretized differently than Eq. [3] Instead of using six neighbors of a lattice 
point in the discretization of the Laplacian, we can use all 26 neighbors of a 3 x 3 x 3 cube around the point. This 
way, we can derive a family of discretizations which is second-order accurate and fourth-order isotropic |24j . The 
Laplacian in this case can be written as 

v^ ijk = (8) 

where D[<pij k \ is given by 

i+i j+i fc+i 

D[(f>ijk] = Y Y ^xyz, (9) 

x— i— 1 y—j~l z=k—l 

and the coefficients Cd only depend on the distance from (x, y, z) to (i, j, k). The values of these coefficients are displayed 
in Table [Tl] for three isotropic discretization schemes [23] . Fig [l] shows a visualization of the distance coefficients on 
three two-dimensional slices in the x direction. The standard anisotropic discretization that we employed in the 
previous section corresponds to c\ = C2 = 0, C3 = 1, C4 = —6. The discrete equation of motion now becomes 



i-1 i i+1 




FIG. 1: The coefficients of the 26 neighbors used in the calculation of the Laplacian at the middle point of the middle two- 
dimensional slice. Table [II] shows the values of these coefficients for three isotropic stencils. 



kjk ~ + ra 2 cj> ijk = 0. (10) 



Again, we assume a solution of the form of Eq. [4] and cancel common terms to get 



2 ; 2 1 2 



(11) 
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coefficient 


Cl 


C2 


C3 


C4 


Scheme A 





1/6 


1/3 


-4 


Scheme B 


1/12 





2/3 


-14/3 


Scheme C 


1/30 


1/10 


7/15 


-64/15 



TABLE II: Three isotropic discretization schemes for the Laplacian. Scheme C uses all 26 neighbors of a 3 x 3 x 3 cube and is 
used in DEFROST. 



with k 2 s given by 



Using e % 



or 

-C2 



c 4 + c 3 [e lk * a + e~ lk * a + e lfe » a + e~ ifc « Q + e lk ' a + e~ lk * a \ 

gi(fc x +fc y )a , ^—i(k x +k y )a . e j(fc;c + fcz)a . e ~i{k x +k z )a _^ e i(k y + k z )a _|_ e —i(k y +k^)a 



_|_gi(fex-fc H ) a _j_ e ~i(k x -ky)a , e i{k x -k z )a __|_ e -i(fc x -fe^)a , e i(k v —kz)a _|_ g — — fca) ra 



+Cl 



2 cos a and cos(a ± /3) = cos a cos /3 =p s i n a sm A Eq. 12 reduces to 



(12) 



;.2 
ft-, 



1 



cir 



(C4 + 2c3 [cos fc x a + cos /c^a + cos fc z a] + 4c2 [cos /c^a cos fc^a + cos fc^a cos k z a + cos k y acos k z a] 



+8ci cos fc^a cos k y acos k z a) 
and the dispersion relation is then 



(13) 



uj (k X7 ky, kz) 



1 



(C4 + 2c3 [cos fc x a + cos k y a + cos fc z a] + 4c2 [cos fc^a cos k y a + cos fc^o cos fc z a + cos k y a cos fc z a] 



+8ci cos k x a cos fc y a cos A; z a) 



(14) 



We perform the same three dimensional simulation as in Section [IT] but now using an isotropic discretization of the 
Laplacian with coefficients c\ ■ ■ ■ 04 given by the discretization scheme C in Table |Ll) This scheme is used in DEFROST 
[19] . We carry out the same initialization and we focus on the same modes as in Section |llj k x = ^(40,28,25) and 
lt2 = ^-(53, 14,2). The results are shown in Table III As in Section |nj the numerically computed mode frequencies 
are different from the continuous ui 2 = k 2 + rn 2 , but agree very well with Eq. 14 In Fig. [2] we plot the numerical 
evolution of (f>(k.2,t) using the anisotropic discretization stencil of Section [TT| and the isotropic discretization stencil C 
in Table |TT| It is clear that the dispersion relation is different for these two different discretization schemes. 



Wavevector 


Lo{h x , ky , hz) 


lo — y/k 2 + m 2 


Numerical 


^(40,28,25) 


4.2139m 


5.4774m 


4.2084m 


^(53,14,2) 


4.0703m 


5.4774m 


4.08m 



TABLE III: Mode frequencies for two modes using the isotropic discretization of Eq. |10| We show the numerically computated 
values, their values in the continuous case ui 2 = k 2 + m 2 and when calculated using Eq. 14 



IV. RADIALLY AVERAGED LATTICE DISPERSION RELATION 



In the previous sections we saw that the dispersion relations on the lattice are not functions of the magnitude 
of the wavevector but rather on its vector decomposition in Fourier space. In some applications however, we are 
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FIG. 2: Numerical evolution of 0(k2,i) with k2 = ^(53,14,2) using anisotropic and isotropic stencils. In the anisotropic 
case w(k2) = 4.2169m and in the isotropic case oj(k2) = 4.08m, which illustrates the mode frequency dependence on the 
discretization stencil used. 



interested in having a radially averaged dispersion relation to match numerical data with theory. One example is the 
computation of the radially averaged two-point correlation function in thermal helds which requires knowledge of the 
radially averaged lattice dispersion relation |21) . In this section we show how to correctly calculate it. 

Given a wavevector magnitude k = |k| we can integrate Eq. 14 over the surface of the positive octet of a sphere 
with radius k, S : k 2 = k 2 + k 2 + k 2 and then divide by the area A = irk 2 /2 to get the radially average dispersion 
relation. First we parametrize the sphere using 



= k 



cos 4> cos 9k x 



k sin <fr cos 9k y + k sin 9k z 



with parameter space Q : < <j> < \tt, < 9 < ~ir. Then the integral of Eq. 10 over the sphere becomes 

LJ 2 (k) 



2 
irk 2 



uj (k x ,k y , k z )da 



2 
7rfc 2 



a/ (k x (</>, 9), k y {4>, 9), k z (</>, 9)) \\N{^ 9)\ \da, 



(15) 



(16) 



where N( 
have \ \N( 



x r' a 



is a normal vector to the surface at the point 



For a sphere of radius k we 



= k 2 cos 9. The expression for the radially summetric dispersion relation then becomes 



! w = -5 - 



,tt/2 ,tt/2 

/ cos 9 / [2c 3 [cos [ak cos <j> cos 9] + cos[afc sin <f> cos 9] + cos[afc sin 9]] 
Jo Jo 



-4c 2 [cos [ak cos <j) cos 9] cos[afc sin <j) cos 9] + cos [ak cos <f> cos 9] cos[ak sin 9] + cos[afc sin <j> cos 9] cos[afc sin 8]] 



+8ci cos [ak cos <fr cos 9} cos [ak sin <f> cos 9} cos [ak sin 9] ] d<j)d9 + 1 



(17) 



Fig. [3] shows the continuous dispersion relation u>i{k) = \Jk? + m 2 and the radially symmetric dispersion relation 
computed from Eq. 17 for two discretization schemes: u>2{k) using the anisotropic stencil with c\ = C2 = 0, C3 = 



1,C4 = —6 and W3(fcjusing the isotropic stencil with c\ = 1/30, C2 = 1/10, c 3 = 7/15, C4 = —64/15. As expected, all 
three agree for low k modes, but they quickly diverge from one another as k gets large. 



V. APPLICATION IN PREHEATING SIMULATIONS 



Even though LATTICEEASY and DEFROST use a discrete grid to perform the fiel d ev olutions, the dispersion 
relation employed is the continuous lo 2 — k 2 + m 2 . As it was shown in Sections In] and III this leads to significant 



G 




FIG. 3: Three dispersion relations: The continuous case with ui\(k) = y/k 2 + m 2 and two radially averaged discretizations, 
u)2(k) and uo3(k). The mode frequency uj and the wavevector magnitude k are given in units of the mass m. 

discrepancies between the lattice frequency and the one predicted in the continuous limit. For this reason, we 
investigate how the use of the lattice dispersion relation affects previous results of preheating simulations. We 
reproduce the results of [5] using the original LATTICEEASY with the continuous dispersion relation and a modified 
version in which it is correctly discretized as shown in Section [IT] We focus our attention on a chaotic inflation model 
with a quartic inflaton potential. The inflaton <p has a four-leg coupling to another scalar field \- The potential for 
this model is 



(18) 



with equations of motion given by 



3-0--V^+(A^+ 5 V 
a a z 



= 



(19) 



ri n ^ 



(20) 



The fields are initialized as Gaussian random fields and the scale factor is evolved self-consistently by the Friedmann 
equations 



8tt 



(21) 



Air 

— (p + 3p)a 



(22) 



where p and p refer to the energy density and pressure of the fields respectively. In LATTICEEASY, the friction terms 
in Eq. [T9| and Eq. [20] are eliminated by appropriate rescalings. For details of the LATTICEEASY implementation we 
refer the reader to the Appendix of [S]. 
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A. Initial Conditions 

The initial conditions are set in Fourier space and then transformed back to get the initial field values in configuration 
space. The simulation starts at the end of inflation and each mode is given a random phase and a gaussian distributed 
amplitude characterized by 



(\fk 



1 

2uj k 



(23) 



where 



w2 = k 2 



'cir 



(24) 



and 



m cff 



d 2 V 
dp 



■9 2 (X 2 



(25) 



for <j> and x respectively. Here () denote spatial averages over the grid. 

Since LATTICEEASY discretizes the Laplacian using the anisotropic scheme of SectiorJIIJ we compare the initial 



power spectrum of the fields using Eq 
of size 128 3 and lattice spacing a 



24 



to what we get using uj 



'cff 



with k 2 s given by EqJ7 We use a box 
0.1 and in order to be consistent with [5] we use A = 9 x 10 -14 and g 2 = 200A. 
All quantities are measured in Planck units (M p = 1.22 x 10 19 GeV) just like in [5]. The initial power spectrum using 
the two different dispersion relations is shown in Fig. [4] for the field </>. Modes with k > 10 show a lack of power 
when using the continuous dispersion relation. This is not surprising if we look at F ig. [3] the radially averaged u>2 is 
smaller than wi for high values of k which results in larger {\<fik\ 
a similar lack of power for high k modes. 



in light of Eq. 23 The power spectrum of x shows 




u? =fe 2 



"ift 



30 

k 



40 



50 



60 



FIG. 4: Power spectrum of field (j> for a chaotic inflation model using LATTICEEASY. The lattice dispersion relation adds 
more power to high k modes. 



Even though the discrepancy in the initial conditions is evident, we must examine whether it has any effect on the 
process of preheating. For this reason, we evolve the fields until the system reaches thcrmalization and compare the 
occupation numbers of the fields for the continuous and discrete dispersion relations. 



B. Occupation Number Density 



by 



The most important variable to calculate during preheating is the comoving number of particles in each held defined 



n f (t) 



(2tt) e 



d 3 kn k {t), 



(26) 



where n k is the comoving occupation number density of particles 



n k (t) 



2uj k 



\fk 



~^\fk\ 2 - 



(27) 



During preheating, the number of particles in each field undergoes exponential growth induced by parametric reso- 
nance. A convenient way to label the end of preheating is by looking at the time when n(t) levels off for each field. 
Fig. [5] shows the evolution of n(t) for both fields </> and x an d matches Fig. 13 of 5J. In order to check how the lattice 
dispersion relation affects preheating, we focus on a time that both fields are long past the exponential growth regime, 
t = 1000 and compare the particle number density. Fig [6] shows the number density for the field <f> at time t — 1000. 
We have performed the run first using the or igin al dispersion relation for both the initial conditions and the 



off 



calculation of the particle number density in Eq. 27 and then using the lattice dispersion relation u> 2 
with fc^ ff given by EqjT] The comparison of the two methods is shown in Fig. [7] for ranges of k in both the low and 
high end of the spectrum which arc most affected by the different dispersion relations. The left part of the figure 
shows that the initial lack of power in high k modes arising from the original dispersion relation has persisted even 
after the end of preheating, giving a consistently wrong number density of the order of 10%. The most interesting 
and unexpected result however is for low k: For the range of 1.5 < k < 3.0, both dispersion relations give almost 
identical mode frequencies, leading to the same initial conditions and definitions for Eq. 27 Yet, even in this case, 
there are features which differ by as much as 10%. This indicates that the lack of power in the initial conditions for 
high k modes affects the process of parametric resonance, leading to an incorrect picture of particle number density 
even long after preheating is over. We therefore conclude that the correct adoption of the lattice dispersion relation in 
preheating simulations is not only desired from a consistency point of view, but required to capture all the available 
features of the underlying theory. 
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FIG. 5: Comoving number of particles as a function of time during reheating for fields <j> and \. Parametric resonance in the 
fields induces exponential growth at early times which then levels off as thermalization is approached. 
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k 

FIG. 6: Particle number density for the field <f> at t = 1000. 




FIG. 7: Lattice dispersion relation effects on the particle number density for the field cj> at t — 1000, shown for small and large 
values of the wavevector magnitude k. 

VI. DISCUSSION AND CONCLUSIONS 

We have calculated the correct lattice dispersion relation for three dimensional simulations employing finite dif- 
ference isotropic and anisotropic discretization methods. We have shown that, on the lattice, the frequency of a 
mode should not be treated as a function of its wavevector magnitude but of its vector decomposition. Moreover, we 
showed how to compute a radially averaged lattice dispersion relation which is important in applications where the 
numerically calculated two-point correlation function is to be matched to the theoretically predicted spectrum. 

Finally, we have shown that the incorrect use of the continuous dispersion relation in numerical simulations of 
preheating after inflation leads to a lack of power in the ultraviolet spectrum of the initial conditions which propagates 
to later times through the evolution of the fields. Consequently, we notice a discrepancy in previous calculations of 
particle number density during thermalization which is accentuated by the incorrect definition of the number density 
operator. 

Even though we have only explored the effects of lattice dispersion relations in preheating simulations, we note 
that modified versions of LATTICEEASY are employed in a variety of other studies including bubble nucleation [25] , 
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gravitational wave production [35] and generation of non-gaussianities [37]. Extra care should be used also in these 
studies to avoid inaccuracies induced by using the continuous dispersion relation on the lattice, particularly in (but 
not limited to) the ultraviolet spectrum. A more detailed study of these effects is left for future work. 
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